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The  Computation  of  Electron  Transfer  Rates:  The  Nonadiabatic 

Instanton  Solution 


Jianshu  Cao,  Camilla  Minichino,*  and  Gregory  A.  Voth 

Department  of  Chemistry,  University  of  Pennsylvania,  Philadelphia,  Pennsylvania  19104-6323 


Abstract 


A  computational  theory  for  determining  electron  transfer  rate  constants 
is  formulated  based  on  an  instanton  expression  for  the  quantum  rate  and  the 
self-consistent  solution  of  the  imaginary  time  nonadiabatic  steepest  descent 
approximation.  The  theory  obtains  the  correct  asymptotic  behavior  for  the 
electron  transfer  rate  constant  in  the  nonadiabatic  and  adiabatic  cases,  and 
it  smoothly  bridges  between  those  two  limits  for  intermediate  couplings.  Fur 
thermore,  no  assumptions  regarding  the  form  of  the  diabatic  potentials  are 
invoked  (e.g.,  harmonic)  and  more  than  two  diabatic  states  can  be  included 
in  the  calculations.  The  method  thereby  holds  considerable  promise  for  com¬ 
puting  electron  transfer  rate  constants  in  realistic  condensed  phase  systems. 
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I.  INTRODUCTION 


Electron  transfer  (ET)  processes  in  chemistry,  physics,  and  biology  have  been  the  subject 
of  a  considerable  number  of  experimental  and  theoretical  studies.  [1,2]  Recent  computational 
approaches  for  computing  ET  rate  constants  range  from  those  based  on  Fermi  s  golden 
rule  [3]  to  explicit  quantum  dynamical  calculations  on  simplified  models  of  ET  processes. 
[4-9]  In  addition,  approaches  derived  from  path  integral  quantum  transition  state  theory 
[10-13]  have  been  developed  to  calculate  the  ET  rate  based  on  the  centroid  density  of  the 
electronic  state  variable.  [14—16]  Despite  the  many  theoretical  and  computational  studies  of 
ET  reactions,  a  unified  computational  approach  has  not  yet  been  developed  which  is  capable 
of  determining  ET  rate  constants  for  arbitrary  values  of  the  electronic  coupling  in  systems 
characterized  by  general  nonlinear  potentials  and/ or  a  significant  degree  of  nuclear  mode 
tunneling.  Significant  progress  towards  this  goal  will  be  described  in  the  present  paper. 

The  underlying  basis  of  the  theory  described  herein  is  the  semiclassical  approximation  to 
the  quantum  partition  function,  [17]  which  can  be  shown  to  be  closely  related  to  the  ther¬ 
mally  averaged  quantum  tunneling  rate  in  metastable  systems.  [18—21]  Along  these  lines. 
Miller  has  suggested  that  the  quantum  reactive  fiux  at  low  temperature  can  be  determined 
by  the  so-called  bounce  trajectory  on  the  inverted  potential  energy  surface,  i.e.,  the  instan- 
ton.  [19]  In  terms  of  the  steepest  descent  approximation,  the  instanton  trajectory  along  the 
periodic  imaginary  time  axis  satisfies  the  Euler- Lagrange  equation,  and  the  quantum  fluctu¬ 
ations  along  the  trajectory  take  the  form  of  a  Gaussian  functional  which  can  be  calculated 
by  evaluating  the  Van  VIeck  determinant.  [22]  The  extension  of  these  ideas  to  the  dissipative 
quantum  tunneling  regime  has  been  discussed  by  Caldeira  and  Leggett  at  some  length.  [23] 
However,  while  the  original  instanton  analysis  is  suitable  when  a  unique  potential  energy 
surface  can  be  assumed,  an  effort  to  include  the  possibility  of  nonadiabatic  transitions  to 
other  potential  surfaces  is  necessary  in  order  to  correctly  describe  electron  transfer  processes 
in  a  general  way.  This  is  the  focus  of  the  present  analysis. 

Many  advances  have  taken  place  in  the  field  of  nonadiabatic  dynamics  simulation  for 
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real  time  quantum  dynamics  (see,  e.g.,  Refs.  [24  27]).  The  theoretical  basis  for  several 
algorithms  [24,27]  is  the  Pechukas  theory  of  nonadiabatic  collisions.  [28,29]  Although  it 
was  originally  formulated  for  real  time  quantum  dynamics,  the  self-consistent  nonadiabatic 
theory  of  Pechukas  bears  a  similarity  to  the  instanton  theory,  both  being  based  on  the 
stationary  phase  (or  steepest  descent)  approximation  to  a  Feynman  path  integral.  [30,31] 
The  former  theory  is  a  real  time  formulation,  while  the  latter  is  in  imaginary  time.  In 
the  present  paper,  the  nonadiabatic  theory  of  Pechukas  will  be  combined  with  the  instanton 
theory  to  yield  a  novel  and  computationally  powerful  approach  for  the  calculation  of  electron 
transfer  rate  constants  under  rather  general  conditions. 

The  present  paper  is  organized  as  follows:  In  Sec.  II,  the  basic  nonadiabatic  instanton 
approach  is  formulated.  A  numerical  algorithm  for  solving  the  equations  is  then  detailed  in 
Sec.  Ill,  and  results  are  presented  for  some  representative  examples.  Concluding  remarks 
are  given  in  Sec.  IV. 


II.  GENERAL  FORMALISM 

To  put  the  formalism  in  the  most  general  context,  we  consider  the  Hamiltonian  for  a 
many-body,  multi-level  system,  given  by 

H  =  H(i{q)  +  Hb{r)  +  Vint{<i,T^)  >  (^•^) 

where  q  is  the  collection  of  N  nuclear  degrees  of  freedom  of  an  electron  transfer  system 
of  interest,  r  is  the  collection  of  the  “bath”  nuclear  degrees  of  freedom,  Hd  is  the  part  of 
Hamiltonian  defined  on  an  electronically  diabatic  basis,  Hb  is  the  bath  Hamiltonian,  and 
Vint  is  the  interaction  potential  between  the  system  and  the  bath.  The  Hamiltonian  Hd  can 
be  explicitly  expressed  in  terms  of  the  elements  ha  (for  the  ith  diabatic  surface)  and  hij  (for 
the  coupling  between  the  ith  and  jth  diabatic  surfaces),  i.e., 

HM)  =  EA.i  +  EE'-v  ■  P.2) 

i  i  j>i 

Here,  the  elements  are  defined  as 
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hn  =  \i)  [K{<l)+Vu{ci)]  {i\ 


(2.3) 


where  K  is  the  kinetic  energy  term  for  the  nuclear  coordinates  q,  and 

hij  =  Vij{q)  (H)(i|  +  li)(il)  , 


(2.4) 


where  the  off-diagonal  coupling  elements  satisfy  the  Hermitian  relation  Vij  =  V*^.  Unlike 
the  usual  adiabatic  barrier  crossing  problem,  the  potential  energy  terms  Vu  in  the  elements 
ha  describe  “simple”  diabatic  surfaces  having,  or  not  having,  potential  wells.  Therefore,  m 
the  most  general  sense  the  quantum  reactive  dynamics  is  induced  by  a  transition  between 
different  diabatic  surfaces  instead  of  taking  place  on  a  single  adiabatic  surface  (e.g.,  a  double 
well  formed  on  the  lowest  adiabatic  potential  surface).  This  formulation  of  the  problem  is 
completely  general. 

Following  a  prescription  originally  proposed  by  Langer  at  zero  temperature  [18]  and  later 
employed  in  various  adiabatic  quantum  rate  calculations,  [20,21]  the  desired  electron  transfer 
rate  constant  can  be  approximated  in  terms  of  equilibrium  quantities  by 


1  ,  InZ 


Tip  Zts 


(2.5) 


with  Zq  being  the  partition  function  of  the  reactant  state,  Z  is  the  overall  system  partition 
function,  and  Zb  is  loosely  defined  here  as  the  “barrier”  contribution  to  the  partition  function. 
The  final  states  are  assumed  to  have  sufl&cient  density  that  ksT  can  be  interpreted  as  the 
rate  of  exponential  tunneling  decay. 

Provided  the  effective  barrier  height  is  significantly  larger  than  the  thermal  energy  in  the 
diabatic  wells,  the  steepest  descent  method  can  be  applied  to  evaluate  the  imaginary  part 
of  the  partition  function  Z  which  leads  to  the  instanton  solution  in  Eq.  (2.5).  A  number  of 
aspects  of  the  instanton  solution  in  various  limits  have  been  elaborated  by  others  (see,  e.g.. 
Refs.  [32-37]).  The  focus  of  the  present  work,  however,  is  on  a  computational  methodology 
to  evaluate  the  instanton  rate  constant  in  the  most  general  case  which  bridges  the  adiabatic 
and  nonadiabatic  (golden  rule)  limits  of  ET.  An  assumption  has  been  made  in  formulating 
this  approach  that  Eq.  (2.5)  is  a  valid  approximation  in  all  limits  of  the  ET  problem.  While 
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numerical  and  analytical  results  presented  below  will  support  this  assumption,  it  has  not 
been  derived  from  first-principles. 

The  stationary  path  of  the  Hamiltonian  in  Eq.  (2.1)  consists  of  the  nuclear  instanton 
trajectory  and  the  self-consistent  electronic  wavefunction  propagation  in  imaginary  time 
arising  from  the  coupling  of  the  two  subsystems.  This  self-consistency  arises  from  the  fact 
that  the  equation  of  motion  for  the  nuclear  coordinates  depends  on  the  imaginary  time 
evolution  of  the  multi-level  wavefunction  which,  in  turn,  depends  on  the  instanton  trajectory. 
A  similar  challenge,  albeit  for  real  time  dynamics,  has  been  encountered  previously  in  the 
study  of  the  dynamics  of  coupled  classical-quantum  systems.  [24,26,27]  Pechukas  was  the 
first  to  provide  a  rigorous  prescription  for  the  self-consistent  stationary  phase  classical-Uke 
trajectory  and  time-dependent  wavefunction  based  on  Feynman’s  path  integral  formulation 
of  quantum  dynamics.  [28,29]  This  elegant  theory  has  since  been  developed  into  various 
approximation  algorithms  for  nonadiabatic  dynamics.  [24-27]  A  Pechukas-like  theory  will 
now  be  developed  for  the  nonadiabatic  quantum  instanton  solution  [cf.  Eqs.  (2.1)  and  (2.5)] 
so  as  to  provide  a  means  for  calculating  the  electron  transfer  rate  constant  under  general 
conditions. 

The  trace  operation  of  the  quantum  Boltzmann  operator  for  the  Hamiltonian  in  Eq.  (2.1) 
involves  a  summation  over  all  the  electronic  diabatic  surfaces  and  an  integration  over  all 
nuclear  coordinates.  Importantly,  however,  this  operation  must  be  rewritten  to  expose  the 
terms  involving  diabatic  state  transitions  which  contribute  to  the  imaginary  part  of  the 
partition  function.  By  inserting  complete  sets  of  diabatic  and  coordinate  basis  states,  the 
tunneling  rate  from  one  diabatic  surface,  denoted  by  [//),  to  another  diabatic  surface,  denoted 
by  \u),  is  related  to  the  following  quantity 

=  J  dr  j  dr'  J  dq  j  dq' 

X  {n,q,r\exp{-pH/2)\u,q^,r'){u,q',r'\exp{-pH/2)\n,q,r)  ,  (2.6) 

where  q  and  q'  are  located  near  the  wells  of  diabatic  surfaces  j^)  and  ]«/),  respectively.  It 
should  be  noted  that  the  two  imaginary  time  propagators  in  Eq.  (2.6)  are  the  same.  (See 


5 


also  Eq.  (C2)  of  Ref  [37],  p.  145.) 

Next,  the  propagator  is  separated  into  the  wavefunction  propagation  of  the  diabatic 
levels  and  the  propagation  arising  from  Hq,  which  is  the  Hamiltonian  excluding  H^,  giving 


=  J  Vr{T)  j  Vq{T)exv{.-So[q{r),T{T)]/ti} 

X  Tf,u[hp,hp/2,q{T)]T^,^[hl3/2,0,q{T)]  .  (2.7) 


Here,  5o[q(T),  r(r)]  is  the  action  functional  excluding  the  contribution  from  Ha,  i.e., 

where  in  is  the  mass  matrix  and  St,  is  the  action  functional  of  the  bath. 

The  quantity  is  the  overlap  between  the  initial  diabatic  state  In)  and  the  final  diabatic 
state  1 1/).  In  an  explicit  form,  the  Bloch  equation  can  be  introduced  to  describe  the  evolution 
of  diabatic  wavefunction,  i.e., 

=  HdhMWr.T')  (2.9) 

SO  that 


t',  q(r)]  =  {u\u{t,t')\h), 


(2.10) 


which  is  a  functional  of  the  system  nuclear  path  q(r)  and  the  imaginary  time  interval  satisfies 

0  <  T  <  h^. 

To  facilitate  further  analysis,  the  bath  average  of  a  functional  /[q(T),  r(r)]  is  introduced 
here  as 

_  /  Pr(rO/[q(r),  r(r)]  exp  {-St,[riT')]/h  -  dT'Vinth{r'),v{T')]/h} 

“  JVv(t')  exp  {-56[r(r0]/a  -  dr'VinMr>)Mr')]/n} 

and  the  quantum  average  over  the  diabatic  basis  for  r  <  7i/?/2  is  denoted  by 

{u\u{hp/2,  T)f{r)u{T,  0)|/j) 


(/(r))d  = 


{v\u{nSl2,T)u{T,  0)|/x) 


(2.12) 


or,  if  r  >  fi/?/2,  then 
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{f{r))d 


(2.13) 


{u\u{hp,  T)f{T)u{T,  hp/2)\fj) 

{u\u{hp,T)u{T,h(3/2)\fi) 

In  Eqs.  (2.12)  and  (2.13),  the  denominators  are  independent  of  the  variable  r  and  /(r)  is 

in  general  a  matrbc.  Both  the  quantum  average  and  the  solvent  average  are  carried  out  by 

assuming  a  particular  nuclear  path  q(T)  and  are  thus  functionals  of  the  nuclear  paths. 

With  the  definition  of  Eq.  (2.11)  in  hand,  one  can  rewrite  the  path  integral  in  Eq.  (2.7) 

as 

Z^,u  =  yX^q('r)exp{-5e//[q(T)]/?i}  ,  (2.14) 

with  the  effective  action  functional  given  by 

SeffHr)]  =  dr  l^q(T)  •m-q(r)  +  W6[q(r)]| 

-n(Zn{r^4;i/?,^/3/2,q(r)]}  +  in{r,.^[W2,0,q(r)]})  ,  (2.15) 


where 


W6[q(T)]  =  -hln^j  rr(r)exp  |-56[r(r)]/n  -  ^  ^  dryint[q(r),  r(r)]/n|  .  (2.16) 

Application  of  the  steepest  descent  approximation  to  Eq.  (2.14)  leads  to  the  equation  of 
motion  for  the  nuclear  coordinates 


m 


^^q(^)  _  /^^d[q(r)]\  ^  IdVintWj)] 


dr'^  \  (9q(r) 


aq(r) 


(2.17) 


which  is  to  be  solved  together  with  Eqs.  (2.9)  and  (2.11)  — (2.13)  to  obtain  the  nonadiabatic 
instanton  solution.  Because  of  the  time  reversal  property  of  the  amplitudes  Tun  T^u^ 
the  instanton  trajectory  is  symmetric  with  respect  to  the  imaginary  time  ^/3/2,  and  so  is 
the  wavefunction.  The  self-consistent  condition  for  the  many-body  nonadiabatic  instanton 
solution  is  two-fold;  the  coupling  between  the  diabatic  states  propagation  and  the  instanton 
trajectory,  and  the  coupling  between  the  bath  distribution  and  the  instanton  trajectory. 

In  order  to  complete  the  instanton  analysis,  the  second  order  functional  derivative  must 
be  evaluated  along  the  instanton  trajectory.  This  procedure  is  numerically  best  implemented 
for  a  discretized  path,  i.e.. 
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(2.18) 


6‘^S  ^  f-tr  ^  -L  \  I  fQ 

_  =  _  (2«„  -  «.,«  -  h,-,)  +  +  'C* 

'a"v;„,iq(T)] 


^3 


+  ^i,3 


dciidcii 


+  e 


where  the  indices  i  and  j  denote  two  different  discretized  imaginary  time  slices,  Qi  and  q^- 
are  the  corresponding  nuclear  coordinates  along  the  instanton  path,  and  e  =  hP/P,  with  P 
being  the  number  of  discretizations.  Here,  C^^ij,  the  bath  fluctuation  correlation  matrix,  is 

given  by 

^  /^VlntlqW]  9Vint[q(T)]\  /  ^Vint[q(^)]  \  /^^mt[q(^)]\  f2  191 

S^K  \  h\  A’ 

and,  Cd,ij,  the  quantum  fluctuation  correlation  matrix,  is  given  by 

_  /9H4q(r)l  ,  .agJq(T)l\  _  /3Hj[q(r)]\  /9gj[q(T)l  \  ^ 

\  dq,  aqs  /,  \  aq,  / j\  aq^  /, 

The  dimensionality  implicit  in  the  above  equations  is  such  that  6‘^S/ dqiSqj  is  a  matrix  of 
dimension  NxP.  When  diagonalizing  this  matrix,  there  will  be  a  negative  eigenvalue  giving 
arise  to  the  imaginary  part  of  the  partition  function,  and  a  zero  eigenvalue  corresponding 
to  the  translationally  invariant  mode.  [20]  The  existence  of  a  zero  eigenvalue  is  an  indica¬ 
tion  of  a  true  instanton  solution.  The  removal  of  the  zero  eigenvalue  requires  the  proper 

normalization,  which  is  explained  in  Appendix  A. 

After  the  preceding  analysis  is  carried  out,  one  arrives  at  the  nonadiabatic  instanton 
approximation  for  the  electron  transfer  rate  constant,  i.e.. 


ksT  — 


(2.21) 


where  W  and  Sinst  are  the  work  and  the  action,  respectively,  along  the  instanton  trajectory, 
and  D  is  a  properly  normalized  determinant  of  the  matrix  in  Eq.  (2.18),  excluding  the  zero 
eigenvalue  (cf.  Appendix  A). 

In  light  of  the  preceding  discussion,  there  are  several  observations  which  can  made: 

(a)  Assuming  a  single  diabatic  surface  in  the  Hamiltonian  [Eq.  (2.1)]  which  contains  a  single 
barrier,  one  recovers  the  well-known  single  surface  instanton  solution.  In  the  case  of  a 
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iiiultil6V6l  systGin,  if  thG  coupling  is  strong  Gnough.  so  thnt  the  nuclGsr  systGin  ulwsys  movGS 
on  thG  lowGst-lying  adiabatic  potential  energy  surface,  the  present  nonadiabatic  instanton 
solution  can  be  shown  to  reduce  to  the  single  surface,  adiabatic  limit. 

(b)  In  the  limit  of  two  weakly  coupled  diabatic  surfaces,  the  Bloch  equation  in  Eq.  (2.9)  can 

be  linearized,  resulting  in  a  transition  amplitude  r',  q(r)]  which  is  proportional  to  the 

off-diagonal  coupling  element  of  the  Ha  matrbe.  This  limit  of  the  theory  thus  recovers  the 
golden  rule  ET  rate  constant.  [3,38] 

(c)  If  the  solvent  is  treated  as  being  classical,  the  bath  paths  r(r)  shrink  to  a  point  and 
Eq.  (2.11)  can  be  rewritten  as  the  configurational  integral 


{f{r))b  = 


/  dr /[q(r),  r]  exp  |-^14(r)  -  dT'Vint[fl{T') ,  r]/^} 


(2.22) 


/rfrexp{-/3n(r)  -  dr’Y^Ur^AIH] 

where  is  the  potential  function  for  the  bath  variables. 

(d)  The  Gaussian  bath  has  a  wide  appeal  in  studying  solvent  effects  in  condensed  media. 
[23,39]  Given  a  harmonic  bath  and  a  bilinear  coupling  between  system  and  bath,  one  can 
explicitly  integrate  out  the  bath  modes  in  Eqs.  (2.11)  and  (2.16),  giving  the  equation  of 
motion  for  the  instanton  trajectory  in  Eq.  (2.17)  as 

where  c(|r  -  r'|)  is  the  imaginary  time  correlation  function  matrix 

nr,  ,  cosh(;i/3a;/2  -  a;|r  -  r'l) 

c(|r-r|)  =  -I  dc-JH - -  ' 

and  J(a;)  is  the  bath  spectral  density  matrix,  related  to  the  elements  of  the  classical  friction 
tensor  riij  (t)  by 

2  /■<»  .  3ij{uj) 


(2.23) 


(2.24) 


naii)  =  -  /  dcj  cos 

^  TT  Jo  u 


uit 


(2.25) 


(e)  In  the  case  of  a  two-state  system  with  a  constant  coupling  between  the  states,  quadratic 
diabatic  surfaces,  and  a  Gaussian  bath,  the  Hamiltonian  becomes  the  spin-boson  model 
which  has  been  often  implemented  in  the  study  of  electron  transfer  (see,  e.g..  Ref.  [39]). 
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III.  RESULTS 


In  this  section,  practical  algorithms  are  described  to  solve  the  equations  in  the  nonadia- 
batic  instanton  theory,  and  numerical  calculations  are  carried  out  for  the  spin-boson  model 
in  order  to  apply  the  theory  to  a  well-known  example.  In  spite  of  its  apparent  simplic¬ 
ity,  the  spin-boson  Hamiltonian  serves  as  the  primary  model  for  investigating  nonadiabatic 
transitions  because  of  its  physical  richness.  Moreover,  the  assumption  of  a  Gaussian  bath 
in  the  spin-boson  model  removes  the  self-consistent  requirement  of  the  instanton  path  and 
the  solvent  distribution,  thus  greatly  simplifying  the  numerical  calculations.  (It  should  be 
noted,  however,  that  there  is  still  the  self-consistent  requirement  of  the  instanton  path  with 
the  nonadiabatic  state  propagation.)  There  is  no  fundamental  problem  associated  with  the 
former  self-consistency  issue  and  a  subsequent  publication  will  deal  with  it  explicitly  for 
multidimensional,  nonlinear  potentials. 

The  major  numerical  effort  in  the  present  theory  is  to  find  the  instanton  trajectory,  that 
is,  to  solve  Eq.  (2.17)  simultaneously  along  with  the  Bloch  equation  in  Eq.  (2.9).  Given  the 
force,  the  equation  of  motion  in  Eq.  (2.17)  is  solved  iteratively  for  a  discretized  instanton 
path.  It  must  be  pointed  out,  however,  that  the  instanton  trajectory  is  neither  a  minimum 
nor  a  maximum  of  the  action,  but  an  extremum  of  the  action.  Consequently,  an  iterative 
method  has  the  possibility  of  converging  the  instanton  in  real  space  to  the  minimum  of  a 
double-well  potential,  which  is  a  trivial  solution  to  the  stationary  condition  in  Eq.  (2.17). 
To  prevent  this  behavior  in  the  iteration  method,  it  is  helpful  to  choose  a  good  initial  input 
trajectory  to  approximate  the  true  instanton  solution.  An  educated  guess  is  the  instanton 
solution  for  the  adiabatic  surface,  which  works  particularly  well  in  the  strong  coupling  region. 
In  the  intermediate  coupling  region,  a  trajectory  solved  for  a  larger  coupling  constant  can 
be  employed  as  an  input  to  the  algorithm.  In  the  weak  coupling  region,  the  adiabatic 
instanton  solution  for  the  cusped  barrier  is  a  good  initial  guess  (cf.  Appendix  B).  The  rate 
of  convergence  depends  on  the  discretization  number  and  the  initial  input.  Generally,  it  has 
been  found  that  around  10^  iterations  will  yield  convergence. 
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Given  a  nuclear  path  q(r),  the  Bloch  equation  Eq.  (2.9)  is  solved  by  numerical  integra¬ 
tion.  At  each  time  step  e  =  the  Hamiltonian  Hd  at  that  time  is  diagonalized  and 

propagated  for  one  step.  The  initial  state  |^)  and  the  final  state  \v)  are  the  right  and  the  left 
diabatic  smfaces,  respectively.  With  the  electronic  wavefunction  in  hand,  one  returns  to  the 
calculation  of  the  instanton  trajectory,  which  in  turn  leads  to  a  new  electronic  wavefunction. 
This  procedure  forms  a  loop  until  self-consistency  is  reached.  In  the  examples  studied  so  far, 
the  convergence  of  the  wavefunction  and  the  nonadiabatic  instanton  trajectory  was  always 
achieved  in  less  than  100  iterations. 

Once  the  instanton  solution  is  found,  the  fluctuation  matrix  of  Eq.  (2.18)  is  computed 
and  diagonalized.  A  vanishingly  small  eigenvalue  will  assure  a  satisfactory  stationarity 
condition  [Eq.  (2.17)]  and  a  negative  eigenvalue  indicates  the  metastability  of  the  particular 
solution  (i.e.,  the  “barrier”  partition  function).  The  prefactor  D  in  Eq.  (2.21)  can  thus 
be  calculated,  and  the  action  S  and  work  W  computed,  hence  yielding  the  instanton  rate 
constant.  In  summary,  the  complete  nonadiabatic  instanton  algorithm  consists  of  following 
steps: 

(1)  An  approximate  instanton  trajectory  is  used  as  an  input. 

(2)  The  stationary  condition  in  Eq.  (2.17)  is  iterated  to  a  converged  trajectory  for  a  given 
electronic  wavefunction. 

(3)  The  Bloch  equation  in  Eq.  (2.9)  is  solved  numerically  for  a  given  nuclear  path. 

(4)  Steps  (2)  and  (3)  are  repeated  until  convergence  is  reached. 

(5)  The  instanton  rate  constant  is  computed  from  Eq.  (2.21). 

As  stated  before,  in  order  to  test  the  method  the  spin-boson  model  was  studied.  In  one 

particular  form,  this  model  is  described  by  the  Hamiltonian 

1  1  N  N 

H  =  -mf  +  A(t^  +  -  a^qof  +  Y.  +  o  Z! 

^  ^  j=i  ^  j=\ 

where  a  is  the  Pauli  spin  matrix,  A  is  one-half  the  tunnel  splitting,  and  the  modes  {a:} 
constitute  the  Gaussian  bath.  The  parameters  were  chosen  in  the  present  case  to  be  ^  =  1.0, 
u  =  1.0,  m  —  1.0,  13  =  5.0,  go  =  5.0.  A  discretization  parameter  of  P  =  200  to  P  =  400  was 


(x^j+uj]x])  ,  (3.1) 
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used  in  the  calculations,  depending  on  the  temperature.  The  parameters  of  the  bath  were 
chosen  so  that  its  spectral  density,  given  in  discrete  form  by  [23] 


^  j=i 


N 


^3 


(3.2) 


reproduced  an  appropriate  friction  kernel  in  the  classical  limit. 

As  a  first  calculation,  a  frictionless  spin-boson  model  was  used  to  verify  that  the  method 
works  in  well-known  limits  and  to  examine  the  numerical  characteristics  of  the  algorithm. 
In  Fig.  1,  the  electron  transfer  rate  constant  is  plotted  as  a  function  of  the  coupling  constant 
A  on  a  logarithmic  scale.  In  the  strong  coupling  region  (i.e.,  large  A)  the  nonadiabatic 
instanton  rate  approaches  the  adiabatic  rate  (dot-dashed  line)  because  the  coupling  is  strong 
enough  that  the  quantum  transition  takes  place  on  the  lower  adiabatic  surface.  In  the 
weak  coupling  region,  the  nonadiabatic  rate  obviously  becomes  proportional  to  the  A^,  as 
predicted  by  the  golden  rule  (dashed  line).  The  golden  rule  rate  in  this  simple  case  is  given 
analytically  by 


ksT 


A^  /7rsinh(6/2) 


tanh(6/4) 

(V4) 


(3.3) 


where  the  activation  energy  is  Ea  =  mu^qQ/2  and  b  =  hPu.  The  adiabatic  tunneling  rate 
reaches  a  plateau,  which  is  the  instanton  rate  for  a  cusped  double- well  discussed  in  Appendix 
B.  It  should  be  noted  that  even  in  this  simple  limit  of  the  spin-boson  model,  the  method  is 
capable  of  dealing  with  an  arbitrary  nonadiabatic  coupling  strength,  bridging  the  adiabatic 
and  nonadiabatic  limits  of  the  ET  dynamics.  It  should  also  be  noted  that  numerically  exact 
methods  exist  for  studying  the  quantum  dynamics  of  spin-boson  model  for  all  values  of  the 
relevant  parameters.  [4-9] 

In  the  adiabatic  limit,  the  instanton  solution  exists  only  in  the  quantum  tunneling  dom¬ 
inated  region,  but  not  in  the  activated  barrier  crossing  region  (for  a  discussion  of  these 
limits,  see  the  review  in  Ref.  [40]).  The  crossover  to  the  instanton  rate  is  given  by  the  well- 
known  criterion  hPui,  >  27r,  with  being  the  adiabatic  barrier  frequency.  A  path  integral 
quantum  transition  state  theory  [10-13]  calculation  can  be  performed  above  the  crossover 
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region  in  the  adiabatic  limit  which  will  bridge  with  the  instanton  solution.  Furthermore,  in 
a  complex  system  all  that  is  required  is  that  a  single  nuclear  mode  be  tunneling  in  order 
for  the  instanton  solution  to  exist.  In  the  nonadiabatic  limit,  the  weak  coupling  induces 
a  nonadiabatic  transition  in  a  small  region  near  the  crossing  of  the  diabatic  surfaces,  thus 
leading  to  a  sharp  barrier  curvature  in  the  adiabatic  surface  which  insures  the  validity  of  the 
instanton  approach.  Therefore,  in  the  golden  rule  region  the  steepest  descent  approximation 
is  always  valid,  even  in  the  classical  limit  of  the  nuclear  coordinates. 

To  further  illustrate  the  characteristics  of  the  nonadiabatic  instanton  solution,  the  follow¬ 
ing  results  are  presented  to  explore  different  aspects  of  the  transition  from  the  nonadiabatic 
limit  to  the  adiabatic  limit: 

(a)  The  nonadiabatic  instanton  trajectories  are  shown  for  A  =  0.01  and  A  =  8.0  in  Fig. 
2.  Obviously,  the  instanton  trajectory  shrinks  as  the  coupling  constant  increases.  On  the 
other  hand,  the  nonadiabatic  trajectory  becomes  independent  of  the  coupling  constant  as 
the  latter  becomes  smaller. 

(b)  Assuming  the  electronic  wavefunction  has  been  determined,  one  can  define  an  effective 
potential  surface  for  the  instanton  trajectory  as 

V„Mr)]  =  (V'l9(T)l)i  ,  (3-4) 

the  derivative  of  which  gives  the  nonadiabatic  instanton  force.  For  comparison,  one  can 
also  evaluate  the  adiabatic  potential  by  diagonalizing  the  Pauli  spin  matrix  in  Eq.  (3.1)  for 
fixed  values  of  the  coordinate  q.  The  effective  potential  is  plotted  along  with  the  adiabatic 
potential  surface  for  A  =  0.01  in  Fig.  3  and  for  A  =  8.0  in  Fig.  4.  As  one  can  observe  from 
Fig.  4,  the  adiabatic  potential  surface  is  a  very  good  approximation  to  the  effective  potential 
surface  for  large  coupling  constant,  whereas  the  cusped  adiabatic  potential  surface  at  the 
small  coupling  constant  in  Fig.  3  is  very  different  from  the  rounded  effective  potential. 

(c)  The  evolution  of  spin  population  is  plotted  for  A  =  0.01  in  Fig.  5  and  for  A  =  8.0  in 
Fig.  6.  As  has  been  stated  earlier,  in  the  adiabatic  region  the  relative  population  on  the  two 
diabatic  surfaces  is  such  that  its  state  vector  forms  the  adiabatic  surface.  In  the  golden  rule 
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region,  the  transition  is  confined  in  a  small  region  near  the  crossing  point  of  the  diabatic 
surfaces  and  happens  rather  dramatically. 

Next,  the  dissipative  quantum  tunneling  region  was  investigated  by  adding  a  Gaussian 
bath  to  the  spin-boson  model  [cf.  Eq.  (3.1)].  The  bath  spectral  density  was  chosen  in  the 
Drude  approximation,  i.e.. 


J{u)  =  r]UJ 


u}t 


up-  +UJ^ 


(3.5) 


where  the  friction  strength  rj  was  1.0  and  the  inverse  of  the  memory  timescale,  Uc,  was  1.0. 
In  Fig.  7,  the  quantum  rate  in  the  dissipative  bath  is  plotted  as  a  function  of  the  coupling 
constant  and  compared  with  the  non-dissipative  rate.  As  is  expected,  the  tunneling  rate  is 
reduced  by  a  substantial  amount  because  of  the  bath  dissipation.  In  addition,  the  dissipative 
suppression  is  stronger  in  the  nonadiabatic  limit  than  in  the  adiabatic  limit. 

Finally,  the  effects  of  anharmonicity  on  the  quantum  rate  constant  were  studied  by 
assuming  diabatic  surfaces  defined  by 


VaiQ)  =  ^"^^^(9  ~  <^zQoT  +  5(9  -  <^290)“*  (3-^) 

where  g=0.01  and  the  other  parameters  are  taken  to  be  the  same  as  in  Eq.  (3.1).  In  Fig. 
8,  the  rate  constant  for  the  frictionless  system  is  plotted  as  a  function  of  the  nonadiabatic 
coupling  constant  A.  Clearly,  introducing  the  anharmonicity  reduces  the  tunneling  rate  and 
the  reduction  is  more  drastic  in  the  adiabatic  region  than  in  the  non-adiabatic  region.  This 
example  illustrates  the  real  strength  of  the  nonadiabatic  instanton  method,  i.e.,  it  is  not 
limited  to  quadratic  diabatic  surfaces. 


IV.  CONCLUSIONS 

In  this  paper,  a  computational  methodology  for  determining  electron  transfer  rates  has 
been  developed.  The  approach  is  based  on  the  instanton  expression  for  quantum  rate  con¬ 
stants  combined  with  a  nonadiabatic  dynamics  formalism  for  treating  the  imaginary  time 
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instanton  dynamics  with  nonadiabatic  transitions.  The  formulation  is  completely  general 
and  thereby  capable  of  treating  nonlinear  diabatic  potential  energy  surfaces  and  multiple 
electronic  states.  It  also  provides  a  computational  method  for  bridging  the  adiabatic  and 
nonadiabatic  limits  of  electron  transfer  processes.  The  theory  was  tested  for  the  well-known 
spin-boson  model,  obtaining  excellent  agreement  with  analytical  predictions  in  both  the 
adiabatic  and  nonadiabatic  (golden  rule)  limits.  In  addition,  it  was  shown  that  both  dissi¬ 
pation  and  nonlinearity  in  the  diabatic  potentials  can  readily  be  included  in  the  calculations 
and  may  have  large  effects  on  the  rate  constant.  In  future  pubhcations,  the  nonadiabatic 
instanton  method  will  be  further  developed  and  applied  to  study  electron  transfer  processes 
in  realistic  systems. 
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APPENDIX  A:  EVALUATION  OF  THE  INSTANTON  PREFACTOR 


In  this  section,  the  prefactor  D  in  Eq.  (2.21)  is  explicitly  expressed  as  a  normalized 
determinant  of  the  matrix  in  Eq.  (2.18).  For  a  free  particle,  the  matrix  describing  the 
quantum  path  fluctuations  is  given  by 


6^S 

SqiSqj 


(Al) 


where  e  =  7i/3/P.  A  normal-mode  transformation  immediately  leads  to  the  eigensolutions  of 
the  matrix  in  Eq.  (Al),  i.e., 


A/  =  2(1  -  cos(27r//P)) 


(A2) 


where  the  index  I  ranges  from  -(P  - 1)/2  to  (P - 1)/2.  Obviously,  I  =  0  gives  a  zero  eigen¬ 
value  which  corresponds  to  the  translational  invariance  of  the  free  particle  space.  Removal 
of  this  zero  eigenvalue  leads  to  the  condition 


HA  = 


(A3) 


which  recovers  the  correct  free  particle  density.  Thereby,  the  instanton  matrix  in  Eq.  (2.18) 
is  normalized  to  the  free  particle  prefactor,  giving 


D  =  lim  -B^det(— -T — 
p->oo  m  oqiOqj 


) 


(A4) 


where  “det'”  stands  for  the  value  of  the  determinant  with  the  zero  eigenvalue  removed.  The 
above  equation  defines  the  prefactor  in  Eq.  (2.21). 
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APPENDIX  B:  THE  INSTANTON  SOLUTION  FOR  THE  CUSPED  POTENTIAL 


In  the  weak  coupling  limit  of  the  spin-boson  model  for  electron  transfer,  the  ground 
state  adiabatic  potential  surface  approaches  a  cusped  parabolic  double  well.  For  such  a 
system,  the  instanton  rate  can  be  exactly  calculated.  [33]  For  simplicity,  a  one-dimensional 
symmetric  double  well  potential  is  considered  here,  given  by 

V{q)  =  ]^rrvJ^{q  -  sign(g)?o)^  (Bl) 


where  the  symbol  ''sign(5)”  stands  for  the  sign  of  q.  The  action  functional  for  a  quantum 
particle  embedded  in  a  Gaussian  bath  then  reads 

5[5(r))  =  dr  +  VI, (r))}  -  L  jf*”  dr'c(|r  -  r>\)  ?(r'),{r)  ,  (B2) 


where  c(|t  -  r'|)  is  the  correlation  function  given  by  Eq.  (2.24). 

The  quantum  rate  problem  for  this  potential  is  most  easily  solved  by  properly  connecting 
the  two  analytical  solutions  of  the  wells  at  the  cusp.  The  resulting  instanton  rate  constant 
is  given  by 


k 


m 


27rnV  ^E(-l)”ar 


:  exp 


mqp 

Yheven 


(B3) 


where  the  factor  is  defined  by 


_ 1 _ 

Ql  +  U!"^  -  f3cn/m' 

Here,  =  27rn/^/3  and  Cn  is  given  by 


(B4) 


Cn  = 

In  the  case  of  a  frictionless  cusped  double-well,  the  rate  constant  k  can  be  expressed  in  a 
closed  form  as 


k  = 


mu  sinh^(6/4) 
h  sinh^/^(6/2) 


exp 


-(3Ea 


tanh(6/4) 

m 


(B6) 
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The  exponential  factor  in  the  above  equation  is  the  same  as  the  one  in  the  golden  rule 
expression  [Eq.  (3.3)],  whereas  the  prefactor  is  by  no  means  the  same.  Thereby,  it  is  necessary 
to  introduce  the  nonadiabatic  coupling  mechanism  in  order  to  obtain  the  correct  limit  for 
the  electron  transfer  dynamics.  The  above  equations,  however,  can  serve  as  a  good  initial 
guess  for  the  nonadiabatic  instanton  algorithm  (cf.  Sec.  3). 
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FIGURES 


FIG.  1.  A  logarithmic  plot  of  the  rate  constant  versus  the  nonadiabatic  coupling  constant  A 
for  the  Hamiltonian  given  in  Eq.  (3.1).  For  comparison,  the  golden  rule  prediction  from  Eq.  (3.3) 
is  plotted  as  a  dashed  line,  and  the  adiabatic  rate  constant  is  plotted  as  a  dot-dashed  line. 

FIG.  2.  The  nonadiabatic  instanton  trajectories  plotted  for  A  =  0.1  and  A  =  8.0  as  a  function 
of  the  imaginary  time. 

FIG.  3.  The  effective  nonadiabatic  potential  defined  in  Eq.  (3.4)  plotted  along  with  the  adia¬ 
batic  potential  surface  for  A  =  0.01. 

FIG.  4.  The  effective  nonadiabatic  potential  defined  in  Eq.  (3.4)  plotted  along  with  the  adia¬ 
batic  potential  surface  for  A  =  8.0. 

FIG.  5.  The  evolution  of  the  population  on  the  two  diabatic  surfaces  plotted  for  A  =  0.01. 

FIG.  6.  The  evolution  of  the  population  on  the  two  diabatic  surfaces  plotted  for  A  =  8.0. 

FIG.  7.  The  dissipative  rate  constant  with  the  bath  spectral  density  given  in  Eq.  (3.5)  plotted 
as  a  function  of  the  nonadiabatic  coupling  constant. 

FIG.  8.  The  rate  constant  in  an  anharmonic  diabatic  potential  given  in  Eq.  (3.6)  plotted  as  a 
function  of  the  nonadiabatic  coupling  constant. 
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